num_cores <- parallel::detectCores() - 1
registerDoParallel(num_cores)Estimating Right-Wing Authoritarianism Values for United States Counties
Using R Regression Modeling
Introduction
The Far Right’s broad resurgence is probably one of the most noteworthy trends in the 21st century’s global political landscape, and additionally politicians representing this shift to the right have gained prominence in not only the US, but in states such as Brazil, India, Russia, Turkey, Hungary, Italy, and Japan. It remains increasingly clear that this trend may be semi-permanent, with over a decade of increasing support for such movements among individuals of all walks of life. “Right wing authoritarianism” is a term not only describing a political identity but a set of beliefs and morals potentially bridging the gap between political science and psychology. A “right wing authoritarianism” scale was developed by psychologist Bob Altemeyer in 1981, describing a set of tendencies that might make one more predispositioned to supporting a right-wing authoritarian leader. These tendencies include high levels of submission to established authorities, willingness to engage in retaliatory actions to protect such authorities, and high degrees of conventionalism and traditionalism, whatever that may be in one’s culture. I will attempt to map such values in the United States, and the county-by-county spatial variance between them. I will take state-by-state or metro-by-metro survey data corresponding with right-wing authoritarian tendencies as defined by Altemeyer, and use a regression model to predict what socioeconomic traits could suggest right-wing authoritarianism in geographic regions.
Perhaps class is a primary predictor of RWA, or urban-rural location, or perhaps Pentecostalism. Modeling RWA by county will allow me to see the factors that correlate the most.
Materials and methods
Firstly, I gathered, clean, and process socioeconomic data (education, income, race, age, and urban/rural) from the census bureau.
I also gathered religion information from the US Religion Census, because I hypothesize that religion could be an additional factor that would explain Right Wing Authoritarianism well.
Then, considering the three major aspects of Right Wing Authoritarianism from Altemeyer, I examined survey data from the PRRI’s American Values Survey and the Pew Research Center. If survey questions closely matched Altemeyer’s definitions, I used them.
It was important that survey data was differentiated on the state level or on some other demographic level (e.g. dividing by income level, race, religion).
I use state-level demographic and religion data as training data, using a generalized linear regression to find the factors that correlate the most with Right-Wing Authoritarianism.
Setting up Parallel Cores for Processing Large Datasets
Loading and Processing Data
I gathered various demographic information from the US Census Bureau API
I began by downloading information on race, education, income, poverty, age, and urban/rural status for all counties in the United States. I used the 2021 American Community survey for all topics except urban/rural status, for which only 2020 Decennial was available. Connecticut recently changed its counties into “planning regions” so there is discontinuity from 2022 onward that will hinder joining to religion data later on. So 2021 has the most recent available data that preserves Connecticut’s old county boundaries.
The Census Bureau API only allows us to download a limited subset of the 2020 Decennial Census online. This does not include urban popuation, for which ACS information is also not available.
So, for information on urban population, I had to download it manually from NHGIS and import it.
The following tables were downloaded:
- B01001: Age by Sex (2017-2021 ACS)
- B01002: Median Age by Sex (2019-2023 ACS)
- B03002: Hispanic or Latino Origin by Race (2019-2023 ACS)
- B15003: Educational Attainment for the Population 25 and Older (2019-2023 ACS)
- B19001: Median Household Income in the Past 12 Months (2019-2023 ACS)
- P2: Urban and Rural (2020 Decennial Census)
The goal is to convert raw values into percentages, and to convert granular data into larger categories I had to convert them all into proportions of the total population. Final column values are to represent: - Total median age - Percentage of total population that is White, Black, Asian, Native American, Pacific Islander, Other Race, Multiracial, and Hispanic - Percentage of those 25 and older who have: received no schooling, been educated up to middle school, high school but no diploma, high school diploma, associate’s degree, bachelor’s degree, post-bachelor’s degree - Percentage of households with incomes: Below 15K, 15-30K, 30-60K, 60-100K, 100-150K, 150-200K, and Over 200K - Percentage of the population in urbanized areas
#Getting a census API key to download census data
census_api_key("65f14507899d7d2b5d92122f92d392799690cb6f", install=TRUE, overwrite=TRUE)
readRenviron("~/.Renviron")#Gathering race data. These are the variables I need to download
acs_racefactors <- c(total_pop="B03002_001",
total_white = "B03002_003", total_black = "B03002_004",
total_indig = "B03002_005", total_asian = "B03002_006",
total_pi = "B03002_007", total_mixed = "B03002_009",
total_hispanic = "B03002_012")
#I convert totals for racial groups in each county into percentages of the population.
acs_race <- get_acs(geography="county", variables=acs_racefactors, year=2021, output="wide") %>%
mutate(white = total_whiteE/total_popE, black=total_blackE/total_popE,
indigenous = total_indigE/total_popE, asian = total_asianE/total_popE,
pacific = total_piE/total_popE, mixed = total_mixedE/total_popE,
hispanic = total_hispanicE/total_popE, total_pop=total_popE) %>%
dplyr::select(GEOID, NAME, white, black, indigenous, asian, pacific, mixed, hispanic)
#The census bureau has information on education level to the year. This is too granular to make an effective model. We need to condense it into broader categories.
educ_factors <- c(total_25 = "B15003_001",
no_school = "B15003_002", pre_k = "B15003_003",
kindergarten = "B15003_004", g1 = "B15003_005",
g2 = "B15003_006", g3 = "B15003_007", g4 = "B15003_008",
g5 = "B15003_009", g6 = "B15003_010", g7 = "B15003_011",
g8 = "B15003_012", g9 = "B15003_013", g10 = "B15003_014",
g11 = "B15003_015", g12 = "B15003_016",
hs_diplom = "B15003_017", ged = "B15003_018",
u_1y = "B15003_019", u_nd = "B15003_020",
assoc = "B15003_021", bachelor = "B15003_022",
master = "B15003_023", prof_sc = "B15003_024",
phd = "B15003_025")
#Mutating to convert into the percentage of the 25-and-over population in seven education cohorts
acs_edu <- get_acs(geography="county", variables=educ_factors, year=2021, output="wide") %>%
mutate(no_school = no_schoolE/total_25E, under_hs = (pre_kE + kindergartenE + g1E + g2E + g3E + g4E + g5E + g6E + g7E + g8E)/total_25E, hs_nodiploma = (g9E + g10E + g11E + g12E)/total_25E, hs = (hs_diplomE + gedE)/total_25E, somecollege = (u_1yE + u_ndE + assocE)/total_25E, bachelor = bachelorE/total_25E, postbachelor = (masterE + prof_scE + phdE)/total_25E) %>%
dplyr::select(GEOID, NAME, no_school, under_hs, hs_nodiploma, hs, somecollege, bachelor, postbachelor)
#Doing the same with income. The census bureau gives 5- to 25- thousand dollar income brackets, we want to convert into broader income categories.
income_factors <- c(households = "B19001_001",
u10k = "B19001_002", u15k = "B19001_003", u20k = "B19001_004",
u25k = "B19001_005", u30k = "B19001_006", u35k = "B19001_007",
u40k = "B19001_008", u45k = "B19001_009", u50k = "B19001_010",
u60k = "B19001_011", u75k = "B19001_012", u100k = "B19001_013",
u125k = "B19001_014", u150k = "B19001_015",
u200k = "B19001_016", o200k = "B19001_017")
acs_income <- get_acs(geography="county", variables=income_factors, year=2021, output="wide") %>%
mutate(u15k = (u10kE + u15kE)/householdsE, u30k = (u20kE + u25kE + u30kE)/householdsE, u60k = (u35kE + u40kE + u45kE + u50kE + u60kE)/householdsE, u100k = (u75kE + u100kE)/householdsE, u150k = (u125kE + u150kE)/householdsE, u200k = (u200kE/householdsE), o200k = o200kE/householdsE) %>%
dplyr::select(GEOID, NAME, u15k, u30k, u60k, u100k, u150k, u200k, o200k)
#The census bureau has age cohort data, divided by sex. Categories are very granular. I think age may be an important variable, but sex, not so much, as the sexes are pretty evenly distributed across US counties. So this is a lot of variables.
age_sex_factors <- c("B01001_003", "B01001_004", "B01001_005", "B01001_006",
"B01001_007", "B01001_008", "B01001_009", "B01001_010",
"B01001_011", "B01001_012", "B01001_013", "B01001_014",
"B01001_015", "B01001_016", "B01001_017", "B01001_018",
"B01001_019", "B01001_020", "B01001_021", "B01001_022",
"B01001_023", "B01001_024", "B01001_025", "B01001_027",
"B01001_028", "B01001_029", "B01001_030", "B01001_031",
"B01001_032", "B01001_033", "B01001_034", "B01001_035",
"B01001_036", "B01001_037", "B01001_038", "B01001_039",
"B01001_040", "B01001_041", "B01001_042", "B01001_043",
"B01001_044", "B01001_045", "B01001_046", "B01001_047",
"B01001_048", "B01001_049", "B01002_001", "B03002_001")
#Aggregating into 5 age cohorts, which represent percentage of the population
acs_age <- get_acs(geography="county", variables=age_sex_factors, year=2021, output="wide") %>%
mutate(age_u18 = (B01001_003E + B01001_004E + B01001_005E + B01001_006E +
B01001_027E + B01001_028E + B01001_029E + B01001_030E)/B03002_001E,
age18_34 = (B01001_007E + B01001_008E + B01001_009E + B01001_010E +
B01001_011E + B01001_012E + B01001_031E + B01001_032E + B01001_033E +
B01001_034E + B01001_035E + B01001_036E)/B03002_001E,
age_35_49 = (B01001_013E + B01001_014E + B01001_015E + B01001_037E +
B01001_038E + B01001_039E)/B03002_001E,
age_50_64 = (B01001_016E + B01001_017E + B01001_018E + B01001_019E +
B01001_040E + B01001_041E + B01001_042E + B01001_043E)/B03002_001E,
age_o65 = (B01001_020E + B01001_021E + B01001_022E + B01001_023E +
B01001_024E + B01001_025E + B01001_044E + B01001_045E + B01001_046E +
B01001_047E + B01001_048E)/B03002_001E,
median_age = B01002_001E, total_pop = B03002_001E
) %>%
dplyr::select(GEOID, NAME, age_u18, age18_34, age_35_49, age_50_64, age_o65, median_age, total_pop)
#Loading in the urban area information separately
urban <- read_csv("data/nhgis_urban_rural.csv") %>%
mutate(urban = U7I002/U7I001) %>%
dplyr::select(urban, GEOCODE)
#Getting geometries for counties
county_geoms <- counties(cb = TRUE, year=2021)
#Joining all the dataframes into one
total_demographic <- left_join(acs_race, acs_edu, by=join_by("GEOID"=="GEOID")) %>%
left_join(., acs_income, by=join_by("GEOID"=="GEOID")) %>%
left_join(., acs_age, by=join_by("GEOID"=="GEOID")) %>%
left_join(., urban, by=join_by("GEOID"=="GEOCODE"))
total_demographic <- left_join(county_geoms, total_demographic, by=join_by('GEOID'=='GEOID'))
total_demographic <- total_demographic %>%
mutate(FIPS = as.numeric(GEOID))The US Religion Census is an organization that tracks religious affiliation across the country. They perform an exhaustive census of every religious organization in the United States to get detailed affiliation data.
There is an API available through OSF in order to download this data, but as of writing the first draft, it appears to be broken. I cannot seem to be able to create an OSF account currently. I do plan on using the API in my final project, but as of now we will just have to make do.
Since the religion census lists affiliation with every single religious organization in the country, the file is very large. I cannot load it directly into the data folder, so I created an abridged version in excel with only the sects I intended on analyzing. These include:
- Total Adherence Rate of All Religious Organizations (TOT)
- Evangelical Protestant (EVAN)
- Black Protestant (BPRT)
- Mainline Protestant (MPRT)
- Catholic (CATH)
- Orthodox Christian (ORTH)
- Mormonism/Latter Day Saints (LDS)
- Judaism (CJUD - Conservative, IJUD - Independent, OJUD - Orthodox, RJUD - Reconstructionist, RFRM - Reform)
- Islam (MSLM)
- Hinduism (HINT)
- Buddhism (MAHBUD - Mahayana, THBUD - Theravada, VAJBUD - Vajarayana)
The Religion Census finds total adherence rates for Christian denominations (e.g. Evangelical Protestant is not a single organization but rather many with similar characteristics)
For smaller religions (Mormonism, Judaism, Islam, Hinduism, Buddhism), I had to manually search through the hundreds of congregations to find these specific groups. There are separate columns for Jewish and Buddhist denominations, but I am combining them, as they are quite small religions.
#Loading the religion data, removing NA values, and creating new columns for Judaism and Buddhism
religion2020 <- read.csv("data/religion2020.csv")
religion2020[is.na(religion2020)] <- 0
religion2020 <- religion2020 |>
mutate(TOTJUDRATE_2020 = CJUDRATE_2020 +
IJUDRATE_2020 + OJUDRATE_2020 +
RJUDRATE_2020 + RFRMRATE_2020,
TOTBUDRATE_2020 = MAHBUDRATE_2020 +
THBUDRATE_2020 + VAJBUDRATE_2020)I also process the religion data by taking NA adherence rate values and converting them to 0. In these cases, it means that there is a negligent amount of congregations of a given sect in a given county, and thus there are not a significant number of adherents.
#Joining US Census Bureau and Religion Census Data
total_demographic <- left_join(total_demographic, religion2020, by=join_by(FIPS==FIPS))So, now I have my demographic baseline data. My next step is to gather survey data from Pew and the PRRI, and predict right wing authoritarianism by county based on these variables.
I created a few maps of different religious affiliations, so we can get a sense of their unique distributions in the US.
It appears that there are many outliers in religious adherence values, with some being suspiciously high, even higher than county populations. This may pose problems for modeling Right Wing Authoritarianism. I will confine values to be within a certain number of z-scores from the mean of all counties.
#Puerto Rico and Guam's county-equivalents are all NA due to the Religion Census not having data for them. So we must remove NA entries.
total_demographic <- na.omit(total_demographic)
#Above three standard deviations is generally considered an outlier. Three standard deviations below the mean for most of these denominations would just result in a negative number. For widespread religions, I will use three. For religions that are found in a limited geographic range, I will use 5
total_demographic <- total_demographic %>%
mutate(EVANRATE_2020 = case_when(
EVANRATE_2020 > (mean(EVANRATE_2020) + sd(EVANRATE_2020)*3) ~ (mean(EVANRATE_2020) + sd(EVANRATE_2020)*3), EVANRATE_2020 < (mean(EVANRATE_2020) + sd(EVANRATE_2020)*3) ~ EVANRATE_2020
), BPRTRATE_2020 = case_when(
BPRTRATE_2020 > (mean(BPRTRATE_2020) + sd(BPRTRATE_2020)*5) ~ (mean(BPRTRATE_2020) + sd(BPRTRATE_2020)*5), BPRTRATE_2020 < (mean(BPRTRATE_2020) + sd(BPRTRATE_2020)*5) ~ BPRTRATE_2020
), MPRTRATE_2020 = case_when(
MPRTRATE_2020 > (mean(MPRTRATE_2020) + sd(MPRTRATE_2020)*3) ~ (mean(MPRTRATE_2020) + sd(MPRTRATE_2020)*3), MPRTRATE_2020 < (mean(MPRTRATE_2020) + sd(MPRTRATE_2020)*3) ~ MPRTRATE_2020
), CATHRATE_2020 = case_when(
CATHRATE_2020 > (mean(CATHRATE_2020) + sd(CATHRATE_2020)*3) ~ (mean(CATHRATE_2020) + sd(CATHRATE_2020)*3), CATHRATE_2020 < (mean(CATHRATE_2020) + sd(CATHRATE_2020)*3) ~ CATHRATE_2020),
TOTRATE_2020 = case_when(
TOTRATE_2020 > (mean(TOTRATE_2020) + sd(TOTRATE_2020)*3) ~ (mean(TOTRATE_2020) + sd(TOTRATE_2020)*3), TOTRATE_2020 < (mean(TOTRATE_2020) + sd(TOTRATE_2020)*3) ~ TOTRATE_2020)
)Loading State- and Metro- Level Socioeconomic Data for Model Training
#Using the same functions used on the county level to extract the same census data for states
acs_race_state <- get_acs(geography="state", variables=acs_racefactors, year=2021, output="wide") %>%
mutate(white = total_whiteE/total_popE, black=total_blackE/total_popE,
indigenous = total_indigE/total_popE, asian = total_asianE/total_popE,
pacific = total_piE/total_popE, mixed = total_mixedE/total_popE,
hispanic = total_hispanicE/total_popE, total_pop=total_popE) %>%
dplyr::select(GEOID, NAME, white, black, indigenous, asian, pacific, mixed, hispanic)
acs_edu_state <- get_acs(geography="state", variables=educ_factors, year=2021, output="wide") %>%
mutate(no_school = no_schoolE/total_25E, under_hs = (pre_kE + kindergartenE + g1E + g2E + g3E + g4E + g5E + g6E + g7E + g8E)/total_25E, hs_nodiploma = (g9E + g10E + g11E + g12E)/total_25E, hs = (hs_diplomE + gedE)/total_25E, somecollege = (u_1yE + u_ndE + assocE)/total_25E, bachelor = bachelorE/total_25E, postbachelor = (masterE + prof_scE + phdE)/total_25E) %>%
dplyr::select(GEOID, NAME, no_school, under_hs, hs_nodiploma, hs, somecollege, bachelor, postbachelor)
acs_income_state <- get_acs(geography="state", variables=income_factors, year=2021, output="wide") %>%
mutate(u15k = (u10kE + u15kE)/householdsE, u30k = (u20kE + u25kE + u30kE)/householdsE, u60k = (u35kE + u40kE + u45kE + u50kE + u60kE)/householdsE, u100k = (u75kE + u100kE)/householdsE, u150k = (u125kE + u150kE)/householdsE, u200k = (u200kE/householdsE), o200k = o200kE/householdsE) %>%
dplyr::select(GEOID, NAME, u15k, u30k, u60k, u100k, u150k, u200k, o200k)
acs_age_state <- get_acs(geography="state", variables=age_sex_factors, year=2021, output="wide") %>%
mutate(age_u18 = (B01001_003E + B01001_004E + B01001_005E + B01001_006E +
B01001_027E + B01001_028E + B01001_029E + B01001_030E)/B03002_001E,
age18_34 = (B01001_007E + B01001_008E + B01001_009E + B01001_010E +
B01001_011E + B01001_012E + B01001_031E + B01001_032E + B01001_033E +
B01001_034E + B01001_035E + B01001_036E)/B03002_001E,
age_35_49 = (B01001_013E + B01001_014E + B01001_015E + B01001_037E +
B01001_038E + B01001_039E)/B03002_001E,
age_50_64 = (B01001_016E + B01001_017E + B01001_018E + B01001_019E +
B01001_040E + B01001_041E + B01001_042E + B01001_043E)/B03002_001E,
age_o65 = (B01001_020E + B01001_021E + B01001_022E + B01001_023E +
B01001_024E + B01001_025E + B01001_044E + B01001_045E + B01001_046E +
B01001_047E + B01001_048E)/B03002_001E,
median_age = B01002_001E, total_pop = B03002_001E
) %>%
dplyr::select(GEOID, NAME, age_u18, age18_34, age_35_49, age_50_64, age_o65, median_age, total_pop)
#Loading a state dataset of urban and rural from NHGIS
urban_state <- read_csv("data/urban_rural_state.csv") %>%
mutate(urban = U7I002/U7I001) %>%
dplyr::select(urban, GEOCODE)
state_geoms <- states(cb = TRUE, year=2021)
#Joining census data
state_demographic <- left_join(acs_race_state, acs_edu_state, by=join_by("GEOID"=="GEOID", "NAME"=="NAME")) %>%
left_join(., acs_income_state, by=join_by("GEOID"=="GEOID", "NAME"=="NAME")) %>%
left_join(., acs_age_state, by=join_by("GEOID"=="GEOID", "NAME"=="NAME")) %>%
left_join(., urban_state, by=join_by("GEOID"=="GEOCODE"))
state_demographic <- state_demographic %>%
mutate(FIPS = as.numeric(GEOID))
#Taking the county-level religion data and using group_by to aggregate it to the state level, weighting by population
state_religion <- total_demographic %>%
group_by(STATE_NAME) %>%
summarize(pop = sum(total_pop),
EVANRATE_2020 = sum(EVANRATE_2020*total_pop)/sum(total_pop),
TOTRATE_2020 = sum(TOTRATE_2020*total_pop)/sum(total_pop),
BPRTRATE_2020 = sum(BPRTRATE_2020*total_pop)/sum(total_pop),
MPRTRATE_2020 = sum(MPRTRATE_2020*total_pop)/sum(total_pop),
CATHRATE_2020 = sum(CATHRATE_2020*total_pop)/sum(total_pop),
ORTHRATE_2020 = sum(ORTHRATE_2020*total_pop)/sum(total_pop),
LDSRATE_2020 = sum(LDSRATE_2020*total_pop)/sum(total_pop),
MSLMRATE_2020 = sum(MSLMRATE_2020*total_pop)/sum(total_pop),
HINTRATE_2020 = sum(HINTRATE_2020*total_pop)/sum(total_pop),
TOTJUDRATE_2020 = sum(TOTJUDRATE_2020*total_pop)/sum(total_pop),
TOTBUDRATE_2020 = sum(TOTBUDRATE_2020*total_pop)/sum(total_pop)) %>%
dplyr::select(STATE_NAME, pop, EVANRATE_2020,TOTRATE_2020, BPRTRATE_2020, MPRTRATE_2020, CATHRATE_2020, ORTHRATE_2020, LDSRATE_2020, MSLMRATE_2020, HINTRATE_2020, TOTJUDRATE_2020, TOTBUDRATE_2020)
#Transforming the religion and census datasets to have the same CRS
state_religion <- st_transform(state_religion, 4269)
state_demographic <- left_join(state_demographic, state_religion, by=join_by(NAME==STATE_NAME))Loading Pew and PRRI Polling Data
Survey Questions Chosen: - PRRI: “Is society becoming too feminine?” (2024) - PRRI: “Is political violence ever justified?” (2024) - Pew: Are there absolute morals? (2024) -
prri_femin <- read_csv("data/prri_femin.csv")
prri_violence <- read_csv("data/prri_violence.csv")
prri_femin <- prri_femin %>%
mutate(FEMIN_MARGIN = FEMIN_AGREE - FEMIN_DISAGREE)
state_demographic <- left_join(state_demographic, prri_femin, by=join_by("NAME" == "State"))
prri_violence <- prri_violence %>%
mutate(VIOLENCE_MARGIN = AGREE_VIOLENCE - DISAGREE_VIOLENCE)
state_demographic <- left_join(state_demographic, prri_violence, by=join_by("NAME" == "State"))
state_dem_na_omit <- na.omit(state_demographic)
set.seed(92322)
fitControl <- trainControl(method="repeatedcv", number=15, repeats=100)
model_femin <- train(FEMIN_MARGIN ~ white + black + indigenous + asian + pacific + mixed + hispanic + no_school + under_hs + hs_nodiploma + hs + somecollege + bachelor + postbachelor + u15k + u30k + u60k + u100k + u150k + u200k + o200k + median_age + urban + EVANRATE_2020 + TOTRATE_2020 + BPRTRATE_2020 + MPRTRATE_2020 + CATHRATE_2020 + ORTHRATE_2020 + LDSRATE_2020, trControl=fitControl, method="rf", data=state_dem_na_omit)
model_violence <- train(VIOLENCE_MARGIN ~ white + black + indigenous + asian + pacific + mixed + hispanic + no_school + under_hs + hs_nodiploma + hs + somecollege + bachelor + postbachelor + u15k + u30k + u60k + u100k + u150k + u200k + o200k + median_age + urban + EVANRATE_2020 + TOTRATE_2020 + BPRTRATE_2020 + MPRTRATE_2020 + CATHRATE_2020 + ORTHRATE_2020 + LDSRATE_2020, trControl=fitControl, method="rf", data=state_dem_na_omit)Variable Importance
#Variable importance for "Society has become too feminine"
varImp(model_femin)rf variable importance
only 20 most important variables shown (out of 30)
Overall
u30k 100.00
postbachelor 69.81
o200k 64.67
u60k 62.39
hs 61.21
asian 58.72
bachelor 54.68
urban 53.07
u200k 38.18
u100k 36.93
EVANRATE_2020 34.04
somecollege 28.35
hispanic 27.46
no_school 27.18
ORTHRATE_2020 24.03
median_age 19.73
CATHRATE_2020 15.24
indigenous 14.01
black 13.92
u150k 13.67#Variable importance for "Political violence can be justified"
varImp(model_violence)rf variable importance
only 20 most important variables shown (out of 30)
Overall
postbachelor 100.00
o200k 96.22
bachelor 92.71
asian 75.41
hs 74.06
u60k 63.33
u15k 51.79
u30k 51.46
urban 50.77
u200k 44.67
no_school 42.25
ORTHRATE_2020 40.14
TOTRATE_2020 29.95
EVANRATE_2020 29.28
u150k 26.87
indigenous 25.27
hs_nodiploma 22.07
LDSRATE_2020 21.76
CATHRATE_2020 19.16
somecollege 18.12Predicting County-Level RWA
femin_predict <- predict(model_femin, total_demographic)
vio_predict <- predict(model_violence, total_demographic)
total_demographic <- cbind(total_demographic, femin_predict)
total_demographic <- cbind(total_demographic, vio_predict)color_fem <- colorNumeric(palette="viridis", domain=total_demographic$femin_predict)
leaflet() %>%
addTiles() %>%
setView(lng=-99, lat=39.5, zoom = 3) %>%
addPolygons(data=total_demographic, color=~color_fem(total_demographic$femin_predict), stroke=FALSE, fillOpacity=0.9, popup=~paste("County:", NAME, "\n Margin:", round(femin_predict, 2))) %>%
addLegend(data=total_demographic, position="bottomright", pal = color_fem, values = ~femin_predict, title="Belief that Society has Become too Feminine")Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'